Metachronal wave and hydrodynamic interaction for deterministic switching rowers 
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We employ a model system, called rowers, as a generic physical framework to define the problem 
of the coordinated motion of cilia (the metachronal wave) as a far from equilibrium process. Rowers 
are active (two-state) oscillators interacting solely through forces of hydrodynamic origin. In this 
work, we consider the case of fully deterministic dynamics, find analytical solutions of the equation of 
motion in the long wavelength (continuum) limit, and investigate numerically the short wavelength 
limit. We prove the existence of metachronal waves below a characteristic wavelength. Such waves 
are unstable and become stable only if the sign of the coupling is reversed. We also find that with 
normal hydrodynamic interaction the metachronal pattern has the form of stable trains of traveling 
wave packets sustained by the onset of anti-coordinated beating of consecutive rowers. 
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I. INTRODUCTION 



Cilia are hair-like extroflections of the cell membrane found in a variety of species from protists to humans, which 
contain active elements (molecular motors, filaments) acting as an internal drive Q. Because of their size and typical 
velocities, the motion of cilia is in the low Reynolds number regime. Ciliary motion can be divided in two stages, 
called power and recovery strokes. The difference between the two is that in the power stroke a higher portion of the 
surface of the cilium pushes against the fluid compared to the recovery stroke so that, as in the breast stroke of human 
swimming, the two effective viscous drags are different, and the filament is able to propel the fluid H,E||i|. Cilia 
normally appear in arrays, and show coordinated wave-like motion, referred to as metachronal wave. The behavior 
of the metachronal wave is thought to be strictly linked to the hydrodynamic interactions between cilia 0, IS El • 
The question of how these collective motions are generated, from the interplay between the internal, active degrees 
of freedom and the external interaction is still open. The scope of this work is to investigate this problem, using a 
simple deterministic model containing very few parameters, consistent with experimental observations and previous 
more detailed modelin g, w hich can be easily simulated and solved analytically in the limit of large wavelengths. 

Modeling of cilia [3, Hfl fill normally requires including the infinite degrees of freedom of an inextensible line-like 
object, its bending elasticity, and its interaction with the fluid (slender body hydrodynamics), plus an active force, 
which can be imposed based on physical observations or treated starting from the "microscopic" internal active degrees 
of freedom [T^.ll3j . The model we present here (section[HJ is extremely simplified and economic in degrees of freedom, 
ft is intended to be treatable analytically. The cilia are represented by point particles, two-state active oscillators 
which we call rowers. The active force is inspired by the switch mechanism introduced by Gueron and others 
Planar or linear arrays of interacting rowers are considered. In a previous work we used the same model to study 
the role of noise |14|. proving that if the switch mechanism of single rowers is purely stochastic the hydrodynamic 
interaction generates metachronal waves which are statistically frustrated by the presence of random fluctuations, but 
can be stabilized by the presence of a short ranged coupling of the internal states, for example of chemical origin. 
An alternative scenario we proposed was that the presence of a coupling between position and transition frequency 
of the single rower would lead to wave like solutions. In this work we would like to pursue this second possibility, in 
the limiting case where the dynamics for the switch is governed deterministically by the configuration of the rower, 
as in the geometric switch of Gueron et alii. 

After an introduction of the model (section [nj , we devote the main body of the paper (section llllj l to the onset 
of metachronal coordination. The discussion is divided in two parts. In the first we discuss an analytical solution 
of the continuum limit of our model equation, which enables to look for the onset of wave-like patterns with large 
wavelengths. In the second part we look at the short wavelengths through numerical simulations. As we will show, 
rowers with a deterministic configurational switch interacting hydrodynamically self-organize in patterns in which 
nearest neighbor particles beat in anti-phase, and propagate trains of wave packets with typical wavelength of a few 
particles. Only with an effectively attractive interaction do long- wavelength wave- like solutions appear. 



II. THE MODEL: ROWERS AND EVOLUTION OF THE INTERNAL STATE 



In the model we adopt, the movement of a single cilium is reduced to that of a low Reynolds number rower that 
maximizes the effective drag in its active phase (power stroke) and minimizes it in the passive, recovery stroke. More 
precisely, a rower is described by two degrees of freedom, of which the first, / is translational and continuous and 
represents a displacement from a reference position. It can be thought of as the displacement of the center of mass 
of the filament from an equilibrium position. The second, a = ±1, is discrete and labels the internal state of the 
object, active or passive. The rowing direction is fixed, while the orientation can in general be left open, to approach 
the problem of symmetry breaking in the generation of fluid flow 0] . The two states carry different effective drags 
-f a = 1 + ecr (with e > 0) to the fluid, corresponding to different surface impacts of the rower (different shapes of 
the filament) in the two phases, together with different potentials (free energy landscapes) V(f,a), that generically 
describe different active or relaxation forces felt by the cilium. This is the implementation of the so-called scallop 
theorem Q at this crude level of description, and makes it possible for a rower to generate a net flow in the fluid. 
There is no interaction between rowers other than the force propagated by the presence of the fluid. This force is 
modeled by the Oseen tensor for low Reynolds number flows, appropriate in the case of cilia [jj. The array of cilia is 
modeled as a linear or planar lattice of rowers labeled by the index n, and the configurations are specified by f n ,o- n . 

This approach to the system contains a radical simplification in the degrees of freedom of the object, a string with 
infinite degrees of freedom, and of the active drive, generated by the collective behavior of many molecular motors. 
This reduction enables us to carry an analytical study. At our level of description, the substitution of cilia with point 
particles does not change qualitatively the interaction induced by the fluid. However, it's a more delicate issue to 
reduce the collective motion of molecular motors to a single dynamic variable. We choose to maintain this variable 



FIG. 1: Potential V(f, a) for k = (+) and k — 1 (o). The dotted lines indicate a = —1, the solid lines a = 1. 



discrete, justified by the experimental observation that the motion of a cilium is divided into two distinct phases and 
by previous, more detailed modeling that suggested this picture [f| 

The evolution of a rower internal state can be modeled generically as a stochastic process 01 i defined by the tran- 
sition frequencies between states. Here we analyze a limiting case where the transition frequencies depend singularly 
on the configuration in a way that reduces the dynamics to a deterministic one. Essentially, a rower contains a switch 
that alters instantaneously its internal state when a particular limit configuration is reached. The dynamics of the 
internal switch is entirely local, in the sense that there is no interaction of chemical origin between nearby rowers. 
With these choices, the evolution of a n may be represented by the following equation containing a Dirac 5 distribution 



dt = -2sgn(/ n (t)) 



df n 



dt 



6(f n (t) - sa n (t)) (1) 



where ±s are the switch points in correspondence of which the discrete internal state is inverted. These parameters set 
the amplitude of oscillation of a single particle and determine its window of motion relatively to the driving potentials 
of the two states (figure HJ. 

Let's now turn to the evolution of the rower displacement and the hydrodynamics. Considering the fact that an 
overdamped motion follows the maximum slope toward the minimum free energy and that we consider no metastability, 
there are generically two possible qualitative choices for the local conformation of the potentials in the two internal 
states. These can be linear, provided the system is far from a minimum, or quadratic if it is close. Therefore, rescaling 
all the constants that are not essential to our discussion (Stokes coefficient, prefactors), we write: 

where the parameter k a £ [0, 1] determines the shape of the potential (figure^. For consistency reasons, here s< 1. 
Thermodynamically, 

- if the switch is close to a minimum of the energy, k a = 1 and we take a quadratic potential, the rower has time 
to dissipate completely its excess energy to the environment before it reaches the switch. The dynamics is a 
cyclical repetition of such relaxation processes 

- if the switch is far from a minimum energy configuration, k a — ► and the potential is linear. The switching 
process is faster than the thermalization time of the rower, that does not have time to dissipate all its energy. 
In this case at the switch the rower must undergo a collision-like process, which conserves the dissipation rate 
and the magnitude of the macroscopic velocity. 

The considerations above refer generically to the active mechanism of model rowers, but leave aside the link with 
real cilia. If one wants to give this drive a microscopic interpretation, it has to be in terms of collective motions of 
the internal motors and elastic degrees of freedom of the cilium. For example, in the linear potential scenario one 
can imagine that motors attach/detach slowly generating a constant force, while in the nonlinear one they attach 
simultaneously, giving the cilium a well-defined minimum energy curvature, and they detach collectively after reaching 
it. In what follows we will set k a = k for simplicity. 

The equation of motion for the rowers has to contain the hydrodynamic interaction. We think of rowers as sources 
for the velocity field and not as boundary conditions, which means introducing a (nondimensional) coupling constant 
a between the fluid and the rowers as a substitute for the geometric constraints, a is proportional to the Reynolds 
number, or to the inverse of the kinematic viscosity. To avoid complications, we do not take into account additional 
boundary conditions, such as walls, but it is straightforward to include them in the model. As the system is in a low 



Reynolds and Strouhal number regime [T^, it is possible to use the regularized Oseen tensor 0] to eliminate the 
velocity field and obtain an evolution equation for the sole rowers degrees of freedom ^ij : 

= (l + ea n (t)) — + 

+ «L n Kra] j- (2) 

where e is the parameter that represents the difference between the effective viscous drags in the two states and 
Q,[n, m] — ^-xx(J + r nm r nm ) is the Oseen tensor projected on the beating direction x of the rowers. Strictly 
speaking, Q[n, m] depends on f n , f m . However, to ease things in an analytical calculation, we approximated it with 
a quantity that depends only on the relative distance of the lattice sites, assuming that the oscillations are small 
compared to this distance. Most of our simulations, though, were carried with the full Q. 



III. METACHRONAL COORDINATION 



The scope of this section is to establish whether hydrodynamic interactions are sufficient to re-phase the rowers in 
absence of a more direct coupling of chemical or mechanical origin. The oscillatory motion of a single rower in an 
array is guaranteed by the structure of its equation of motion. A mean field description of the array can be carried out 
considering the overall effect of the velocity field generated on one rower by all the others. This procedure is outlined 
in appendix ^ and leads to the main result that a collection of rowers can generate a macroscopic flow if e ^ 0, and 
that, provided there is no intrinsic orientation in the beating mechanism, symmetry will be broken. However, in a 
description that goes beyond mean field, nothing can be said about the beating time, which can in general vary with 
the dynamics, so that the question can be restated as whether this variability in the beating time is stabilized or 
disrupted by the hydrodynamic interactions. The problem is hard to approach analytically due to the discontinuous 
nature of the switch a and the nonlincarities. However, the continuous limit of the model, which describes the long 
wavelength behavior of the system, can be approached analytically. The results obtained in this way can then be 
compared with the numerical study of the discrete case. 

From the point of view of solving the equation of motion one has to investigate: 

(1) The existence of wave-like solutions. We define metachronal waves motions of the kind f n (t) = f(t ± rx n ), and 
simple metachronal waves those for which / (and thus a) is a periodic function. 

(2) The stability and attractivity properties of these solutions. 

(3) Their statistical weight in a macroscopic description of the system on large time scales. As one cannot establish 
a priori its initial conditions, the metachronal solutions will be significant if the phase space volume of the 
initial conditions they attract is nonzero, and the relaxation time scales do not exceed a cutoff defined by the 
lifetime of the system. 

In what follows we will be mainly concerned with the first two points. The third point will be approached in general 
numerically for systems of a few rowers (short wavelengths). 



A. Metachronal pattern in the large wavelength limit. Continuous model 



We will show that in the continuum limit of the equation of motion [2] it is possible to find simple metachronal wave 
solutions and study their stability analytically. We can take the continuous limit analyzing selectively metachronal 
solutions whose wavelength is large compared to the spacing between rowers. Then f n (t) = /(x n ,t) becomes the 
continuous field /(x, t) and we can rewrite the hydrodynamic interaction tensor fi[n, m] as (—V 2 + q 2 )^ 1 , where we 
incorporate also the possibility of screening with inverse screening length q. With one inversion of this operator, the 
evolution equation for the continuous / can be written as 



(<z 2 - V 2 



^ + (l + ea)(kf-*) 



= -a(kf - a) 



(3) 



The laplacian in the above expression is one dimensional, along the fixed direction of beating x . We look for planar 
metachronal wave solutions with the ansatz / = f(t — tx), so that we can reduce to a 1+1 dimensional problem. 



The transverse hydrodynamic interactions are irrelevant as the rowers are constrained to beat in one dimension, and 
the anisotropic terms can be absorbed in the prefactor of the interaction tensor. Calling y — t — tx, we can assume 
without loss of generality that y = is the coordinate of a wave-front, where the switch a has a jump. This translates 
in the condition 

a = 0(-y) - 9(y) 

where 9 is the Heaviside step function. The local displacement / can be decomposed in the sum of power stroke and 
recovery stroke parts, f+0(—y) + f-8(y), this implies that 

vf = f+6(-y) - /-%) 

With this procedure one obtains two linear third order differential equations for /_ and / + , together with four 
joining conditions for f± and their derivatives in correspondence with the jumps of a. Metachronal solutions can 
be constructed starting from the initial condition = and generating a succession of wave-front coordinates 
yW, j/ 2 ) . . . z/™) imposing the joining conditions above on the solutions of the third order differential equations for 
f±. The iteration of this process, starting from the initial conditions for f± and its first and second derivatives, 
or equivalcntly on the vector of the (two at the most) independent arbitrary constants (A±,B±) of the solution of 
the differential equations, generates a flux in phase space, described by an affinc transformation. The existence of 
a fixed point of the succession (A^ , B±^ ) and its attractive properties determine the nature and stability of the 
constructed metachronal wave. At every iteration, the y^ must satisfy the relation /(j/™-*) = ±s. Note that despite 
every step involves linear operations, a strong nonlinearity is introduced by the inversion for i/™) of the solutions 
of the differential equations. It is also important to stress that the solutions constructed in this way are in general 
not periodic, and may have a domain of existence which is bounded in y, as after a number of iterations it could be 
impossible to invert for ?/ n ' . More details on these calculations for the exemplifying case k — 1, e = are reported in 
appendix [5] 




FIG. 2: Analytical solution of the continuum equation computed for the case s = 0.8, p = 0.01, k = 1. The fixed point 
values of the parameters in this case are A± = —0.435, B± — —1.456. 

Our results can be summarized as follows. We find that, from the point of view of the existence and stability of 
metachronal solutions, e plays no qualitative role, so that the problems of flow-generation and syncronization can be 
separated. For this reason, the discussion is independent on e and can be simplified restricting to the case e = 0. A 
necessary condition to find that s < 1, which gives a good consistency test. 

- For quadratic potentials (k > 0) the qualitatively relevant parameter is p = ^r-. Fixing s, there is a critical 
value p c such that for p > p c there exist no solution. This set an upper critical wavelength A c for the metachronal 
waves. For < p < p c , a fixed point exists, and the stability analysis gives a marginally stable saddle point 
in the planes of coefficients (A±,B±). In other words there is one line (a region of zero measure) of stability 
in this plane, with runaway hyperbolic trajectories around it. We report this case in figure [3] as an example. 
Finally, for p < the wave like solutions are always stable. However, this last condition implies p < 0, therefore 
an effective "attractivity" of the hydrodynamic interactions. We will discuss below two cases that could lead to 
this situation. 

- In the case k = (linear potentials) the solution always exists. Stability analysis gives solutions that are always 
unstable (two eigenvalues > 1) for p > and stable for p < 0. 

In both cases we can conclude that metachronal solutions at long wavelengths are unstable unless p < 0, and the 
interaction becomes attractive. We would like to spend a few words on the possible physical meaning of this change 
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FIG. 3: Stability analysis for s= 0.3, p = 0.05, k = 1. The fixed point values of the parameters are A± — ±0.337, B± — 
±0.865. The fixed point is marked with a black square, and the lines with arrows indicate the direction of evolution of 
different initial conditions, showing that the fixed point is marginally stable. The trajectories are obtained iterating the 
transformation on (A±,B±) starting from different initial conditions. The lines reported are Bezier interpolations of these 
(discrete) transformations. 

in sign, with an effectively attractive interaction. One first consideration is that interaction between colloidal objects 
can be more complicated than how we represent it, in presence of hydrodynamic effects or charge. For example, it 
has been speculated that the presence of a wall in combination with a surface charge, all neglected in our model, 
lead to effectively attractive potentials between colloids of like charges . Lubrication forces are another possibility, 
provided that cilia are able to get sufficiently close to each other. For rowers, this last condition is contrary to the 
assumption of small oscillations in the theoretical calculation but can be tested through simulations with the full 
Oseen tensor. Coming back to the more orthodox case p > 0, the results above show that in the continuous limit 
simple metachronal waves (with nearest neighbor rowers in phase) exist for all wavelengths (below a characteristic 
one, for k ^ 0), but they are (marginally for k ^ 0) unstable. These results are not sufficient to establish in general 
the statistical behavior of metachronal waves, as the hypotheses adopted here reduce the analysis to a particular class 
of solutions and the local region of the phase space that surrounds them. For negative p only we can conclude that 
wave-like solutions attract all the initial conditions. In the other cases, where the solutions are not attractive, we have 
to look for other basins of attractions in the phase space. One possibility is that these lie in the shorter wavelength 
region, which is overlooked by the continuous model. This motivates the numerical analysis of the next sections. 

B. Numerical simulations with many rowers 

A confirmation of these results, and more insight on the behavior of the system comes from numerical simulations 
of linear and two dimensional arrays of many (50-250) rowers. These were run starting from random initial conditions, 
with nearest neighbor Oseen-like interactions, the simplified interaction tensor £l[n, m] and the full f2[/ n , f m ], robustly 
showing the same phenomenology in all cases. In particular, for linear arrays with k > 0: 

- For p > the relaxed solutions look like trains of traveling wave packets where nearby rowers show anti- 
coordinated motion. In the patterns observed in the simulations the packets have typically a characteristic 
length of the order of 10 particles. Their traveling direction coincides almost always with the beating direction 
of single rowers. 

- For p < long wavelength traveling solutions appear. These always have a wavelength that is exactly the size 
of the array. For every solution of this kind traveling in one direction, there is another one in the opposite one, 
leading to mirror symmetric standing waves. The same solutions are found if a is kept positive, and a nearest 
neighbor attractive interaction is added. 

This is exemplified in figure^] For linear arrays with k = wave-like solutions appear only for p < 0. Two dimensional 
arrays show the same qualitative behavior. This makes the wave patterns propagate in directions that are different 
from the beating direction of the rowers (non-simplcctic, in the language of metachronal waves). 

We can give a heuristic argument based on the switch mechanism to account for this anti-phase coordination. Let's 
consider three consecutive rowers X,Y,Z, and suppose they are in phase. Immediately after X reaches the ±1 switch it 
feels a strong negative force due to the change in the potential, which is propagated mainly to Y and much less to Z. 



(a) (b) (c) 

FIG. 4: Metachronal solutions from numerical simulations of linear arrays of rowers interacting with the full Oseen tensor. The 
distance between successive rowers is 2 in our rescaled units and e = 0.25. (a) and (b): case p > 0. a = 0.4, k — 1, 150 rowers. 
In (a) a part the configuration /„ as a function of the rower index is shown. The dashed line shows the actual configuration, 
while the solid one, obtained connecting the points relative to odd rowers, outlines the shape of the wave packet, (b) represents 
the time evolution /(t) of the 40th rower of this array, (c): case p < 0. a — — 0.1, k = 1,200 rowers. The configuration f n 
shows a long wavelength metachronal wave. The part of the configuration for i > 100 is mirror symmetric, due to the same 
pattern traveling in the opposite direction 

Provided Y was going toward the same switch, it will be slowed down, and therefore de-phased with respect to X. If we 
imagine now that X,Y,Z are in anti-phase and A reaches the +1 switch. Now the force felt by Y (and much less by Z) 
will be driving it toward the -1 switch, reinforcing the dephasing. As a result, the anti-phase beating between X and 
Y is more stable than the coordinated motion. According to this argument, this will be the case in all the instances 
where the influence of nearest neighbors is stronger than that of far away ones. More formally, one can find a reason 
for this behavior in the analytical solution of the continuum equation. Configurations where rowers, interacting with 
the normal Oseen tensor have opposite phase are overlooked by the continuum limit. However, roughly speaking, 
these configurations lead effectively to a change of sign of a in the continuum limit equation, as one can easily see 
imposing / ovo „ = — / odd in equation[2| Therefore one can think that the stable waves found for p < in the continuum 
model have a trace of this anti-phase behavior. 

C. Stability of the metachronal wave at short wavelengths. 

To analyze the system at shorter wavelengths we ran numerical simulations of linear arrays of a few rowers with 
periodic boundary conditions, choosing to truncate the interaction to nearest neighbors. We interpret the boundary 
condition as a constraint on the wavelength of the resulting wave patterns. For up to four rowers, we analyzed the 
system using Poincare maps. This study leads to more general considerations on the statistical weight of metachronal 
solutions. The maps are obtained graphing the positions of two rowers when a reference one reaches the switch over 
many cycles of the motion and for different initial conditions. If the resulting graph is, or converges to a single 
point, the motion is coordinated. If it is a closed orbit the coordination is quasi-periodic, which means that the 
motion in general does not look like a traveling wave. Finally, if the resulting graph is a random scatter plot, the 
motion is chaotic. Indeed, the system can be either chaotic or quasi-periodic if the potentials are linear, (k — 0), 
as can be seen in figure |SJ The phase space volume of the quasi-periodic region delimited by the concentric loops 
is therefore a measure of the statistical weight of quasi-periodic solutions. The volume of this region decreases with 
increasing a. The case k ^ is somewhat different. Here, there is no chaos, and the stable tori become attractors 
to a fixed point of the parameter space (figure El (a)). This means in principle that all the initial conditions lead 
to pattern formation. Therefore, the relevant parameter becomes the relaxation time. This diverges as k — > as a 
power-law with exponent b ~ 0.85 and as a — > with exponent c ~ 0.74. The curvature of the potential, which as we 
discussed has to do with the activity of the microscopic active degrees of freedom, seems therefore to be important in 
determining the organization properties of the system. On the other hand, the phenomenology observed for different 
values of a is qualitatively consistent with what observed experimentally for arrays of cilia beating in fluids with 
varying viscosity 0, El • 
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FIG. 5: Poincare map for the case k — 0, three rowers, periodic boundary conditions. Here, e = 0.02 and s= 0.9. The allowed 
configurations are confined in the square region [-s,s] x [-s,s]. However, to increase clarity, the figures were divided in four regions, 
corresponding to the different values of the two switches, which were "folded outside" , so that the axes correspond to the actual 
configuration of the rowers after a mirror reflection dependent on the state of the switches, (a) a — 0.1. (b) a = 0.01. 
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FIG. 6: (a) Poincare maps for k 7^ 0, three rowers, e = 0.02. Left: a = 0.01 attracting trajectories for different values of k. A 
k = 0.15. Ok — 0.2. Ok — 0.3. O k = 0.4. Right: different initial conditions for k = 1. (b) log-log plot of the relaxation time 
as a function of k for a — 0.01(D) and a = 0.1 (O). The power law fits (solid lines) yield an exponent of 0.85. (c) Relaxation 
time as a function of a for k — 0.3. The power law fit gives the exponent 0.74. 

IV. OVERVIEW AND CONCLUSIONS 

We have presented a simple model system of two-state low Reynolds number oscillators called rowers as a generic 
framework for the problem of cooperation of cilia. The dynamics adopted in this work, specified setting the transition 
rates between the two potentials, is entirely deterministic, determined by a switch mechanism coupled to the configu- 
ration. We solved analytically for wave like solutions the continuum, long wavelength, limit of the equation of motion 
for an array of rowers with hydrodynamic interaction and we analyzed the stability of the solutions, confronting with 
results from numerical simulations. Finally, we analyzed through Poincare maps the phase space dynamics of systems 
of a few rowers, to study their behavior at short wavelengths. 

Our most important result is that metachronal patterns exist at all wavelengths (below a characteristic one, for 
k ^ 0), but long wavelength solutions are (marginally for k 7^ 0) unstable. The stable patterns have the form of 
consecutive wave-packets where nearest neighbor oscillators are in anti-phase, propagated with constant speed, with 
a characteristic length of a few rowers. We showed that the statistical weight of these solutions can be determined 
numerically imposing an upper cutoff on the wavelength of the pattern. Only in the presence of a reversed coupling 
constant, can long wavelength metachronal solutions be stable. We proposed two possible physical reasons for this 
reversal in sign. Deterministic switching rowers, as two state oscillators, show a rich and unusual phenomenology, 
of which we could explore a number of aspects. Their behavior is in many ways opposite to our usual notion of 



oscillations, starting from the fact that no normal modes can be defined, but the oscillators self-tune to a chosen 
frequency determined by the characteristic relaxation times in the two states, much as in systems close to a Hopf 
bifurcation [lflj . 

Comparing the behavior our abstract entities with that of real or model cilia, the first puzzling question seems to 
lie in the anti-phase motion. As discussed, a solution of this could lie in a short ranged interaction with a different 
origin. One good candidate for this are lubrication forces, as real cilia can be really close to each other. Also, a short 
ranged synchronization between the switches of chemical origin could lead to the same result, consistently with the 
scenario proposed in our previous work. The relevant parameters in our discussion are the stiffness of the potential k 
and the hydrodynamic interaction coupling strength a. The first is related to the internal active degrees of freedom, 
which are hard to access experimentally, while the second can be used for a qualitative comparison of our results 
with experiments where arrays of cilia are observed beating in fluids with varying viscosity 0, Q • One other question 
is the relation with more detailed models of cilia and their internal drive, in particular with the geometric switch 
model of Gueron and collaborators [.5, 6J. Rowers, with their few degrees of freedom constitute a system much more 
under control than filaments to test. We can conclude that generically simple hydrodynamic interaction does not 
synchronize but anti-synchronizes nearest neighbor rowers, so that, if filamentous objects are to be synchronized by 
a similar mechanism, an extra (to be found) ingredient is needed. 
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APPENDIX A: MEAN FIELD APPROACH, SYMMETRY BREAKING 



It is possible to estimate the magnitude of the characteristic beating times and the macroscopic speed generated on 
the fluid without having to solve explicitly the equations of motion, using simple self-consistency considerations |l4j . 
This is done writing the equations of motion QandEl for a single rower in a constant external effective velocity field 
v along its beating direction and supposing that this velocity is generated by the effect of the surrounding rowers. 
Equation [3 looks then like 



dfojt) 
dt 



d[-V(f Q (t),a (t))} 
dfo 



where we labeled conventionally with the index the rower around which we do the selfconsistent calculation. It is 
then straightforward to calculate the beating times t± for this single rower, defined as the times required to go from 
-s to +S and back respectively: 



U 



(l+e)(l+kS)+v 

(l + £)(l-fcS) + K 



t- 



[k(l - i)\- y ln 



(l- £ )(l-fcS)- 



(l-e)(l+fcS)-o 



These quantities can be used to determine self-consistently the absolute value of the "macroscopic" fluid velocity v, 
taking into account the average force exerted by the single rower on the fluid in one cycle. 



Vsc 



dV(f,a) 



1 



T Jo df 



dt 



where the total period T = t + + i_ depends on v. Therefore 



^sc = ~™ 



n/0 



0,n 



VSC 



t- 



1 + 6 



4es 



(Al) 



In this expression the summed hydrodynamic propagator simply plays the role of a multiplicative constant, and s 
determines the self-consistent value for the velocity (see figure 0). For e = 0, v$c = 0, and the scallop theorem is 
found again in the description of a collection of rowers. The same argument yielding a non zero self-consistent velocity 




FIG. 7: Selfconsistent velocity calculation. The right hand side (solid line) and the left hand side (dashed line) of equation 
IA1I are plotted as a function of v in the graph above. The intersection between the two lines yields vsc- The values of the 
parameters for the curves shown are k = 0.8, e = 0.25, s = 0.8, and a X^ n? ^o ^°, n = 0.2. With these values, wsc — 0.365 in 
nondimensional units. 



can be extended to the case of one - or a collection of- rowers whose rowing direction is not fixed. This is obtained 
taking reflection symmetric potentials V(f, a) — V(—f,a), one of which with a single minimum for / — 0, the other 
with a double well. A thermal noise needs to be added for the problem to be well posed. Rowers are able to break 
spontaneously the symmetry to generate a flow in the fluid 14]. For real cilia, the problem of symmetry breaking can 
be relevant in the context of generation of left-right asymmetry through nodal flow in vertebrate embryos |17| . 



APPENDIX B: EXAMPLE OF EXPLICIT SOLUTION OF THE CONTINUOUS MODEL AND ITS 

STABILITY 

We will solve equation |3 analytically with the ansatz / = f(x — rt) on the solution. For simplicity we can restrict 
ourselves to the case k = 1, q = 0, e = 0, as the general case carries no further conceptual complication. Calling 
y = x — rt, the equation reads: 

- T 2 (f" + /") + a(f tV = (Bl) 

where ' indicates derivatives with respect to y. For a transition of a from 1 to -1 at the wave-front y = 0, the right 
joining conditions are, as already discussed, 

o = 9(-y) - 9(y) (B2) 
f = f + 0(-y)+f-6(y) (B3) 
erf = f+0(-v) ~ f-0(y) (B4) 

where 9 is the Heaviside step function, and f± = f±(y)- Analogue expressions hold for the transition — 1 — ► 1. The 
decompositions above generate two linear ordinary differential equations for f±. 

f±+f±-pf± = TP (B5) 



with p — Moreover, the same conditions IB2I IB3I IB4I and their derivatives can be substituted in equation IB II 
obtaining an expression containing terms in 9(±y), 8(y) and its derivatives. Equating all the terms to zero one obtains 
three joining conditions. That is, 

(/+-/-)lswitch = °; ( B6 ) 

(/+-/-) I switch = 3! ( B? ) 
4 

(/+-/-) I switch = -g ( B8 ) 



Here, and more generally in the case e = 0, the two equations IB5I are the same with the identification — /_ = /+. 
The solution of equation IB5I is easily obtained; f± = ±1 is always a particular solution, and one has to solve the 
characteristic equation z 3 + z 2 = p. This admits three real solutions for < p < p, and one for p < O.p > p, p = 4/27. 
We will analyze in detail here the case with three solutions {—Zi, —Z2, 23), with Z\, z%, Z3 > . In this situation, 



/±(y) = ± (1 + A±e- ZiV + B±e- Z *y + C ± e Z3y ) 



The constant C± can be eliminated using the condition /±(0) = =ps, meaning that after the first jump the rower is 
located at the switch. The next step is to evolve this solution up to a certain y± where the next switching event will 
take place, imposing that 

f±(y±) = ±b (B9) 

y± is obtained inverting this last expression, and has to satisfy the joining conditions IB8I for the next "piece". For 
example, supposing we start from state a = 1 

(f' + ( y+ ) f'_(0)) = \ ; (f'i(y+) £(0)) = ~ 

This gives a linear transformation (A + ,B + ) — > (A-,B-). A complete solution can be constructed iterating this 
procedure. This solution is in general non-periodic, as y± may vary at every step, and also cease to exist. The 
equations that determine y+ and yi^ at the n-th step can be written as 

4 n) #(i,3)(2/i n) ) + 4 n) #(2,3)(y£ n) ) - (1 + s)E 3 (y^) + (1 - s) = (BfO) 
where we used the notation 

H (ii3) (y) = e-™ - e-™ i = 1, 2 



Ei(y) = e-« v i = l,2; E 3 (y) = e 

The joining conditions for step n are 



A^=E 1 (yp ) )A^+w 1 
4 n+1) =E 2 (y^)B^+w 2 



where w±,W2 are rational functions of the solutions zi, ziz 3 . 

A simple (periodic) metachronal solution exists when the transformation has a fixed point, which can be imposed 

setting the equality between A^ +1 ^ and This wave is characterized by a unique y+ = y_ = y solution of 

Trm*™® + Trm*™® ~ (1 + s)Ea{g) + (1 - s) = 

and by coefficients A_ = A + = A, B- = B + = B given by 

A= 

1 + e^ z iy 



l + e -z 2 y 

The equation for y admits solution only if s < 1 , and, for any given value of s, if p is lower than the critical value p c < p 
introduced in the paper, which can be found numerically. This sets a maximal wavelength A c for the metachronal 
wave. The stability of the solution can be evaluated linearizing the flow, starting from the point (^4 + 6 A, B + SB) 
in parameter space, inverting equation IB 1 01 for y and calculating the total variation (5A,5B) from IB ill In this case 
this yields one negative and one positive eigenvalue corresponding to a marginally stable fixed point. The procedure 
outlined in this appendix can be carried out more in general, leading to the results discussed in the body of the paper. 
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